library("copula", lib.loc="~/R/win-library/3.1")


setClass("AEP.frame", representation(
  data = "data.frame",  
  nb_obs = "numeric",
  nb_var = "numeric",  
  N = "numeric",  
  binary_matrix = "array",
  sum_ones = "array",  
  alpha = "numeric",
  iteration = "numeric",
  m = "array",
  bandwidth = "numeric",
  data.weighted_bandwidth = "matrix"
  )
)

setMethod("initialize", "AEP.frame", function(.Object, data, alpha=2/(dim(data)[2]+1), iteration=1, bandwidth=rep(0, c(dim(data)[2]))) {
  
  if (missing(data)) stop("syntax error")
  stopifnot(is.data.frame(data), is.numeric(alpha), is.numeric(iteration), is.numeric(bandwidth))
  
  if (dim(data)[2] == 1) stop("'data' must have at least two columns" )
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")   		  
  if (iteration < 1) stop("the parameter 'iteration' must be >= 1")			  
  if (length(bandwidth) != dim(data)[2]) stop("the parameter 'bandwidth' must be of the same length as the number of column of 'data'")
  
  .Object@data <- data
  .Object@nb_obs <- dim(data)[1]
  .Object@nb_var <- dim(data)[2]  
  .Object@alpha <- alpha	
  .Object@iteration <- iteration
  .Object@N <- 2^.Object@nb_var
  .Object@binary_matrix <- array(0, c(.Object@N, .Object@nb_var))			
  .Object@sum_ones <- array(0, c(.Object@N))  
  .Object@bandwidth <- bandwidth  
  .Object@data.weighted_bandwidth <- t(.Object@bandwidth * t(.Object@data))  
    
  #Initialize the matrix of binary and the sum of ones in each line
  for (i in 2:.Object@N) {
    .Object@binary_matrix[i,] <- .Object@binary_matrix[i-1,]
    j <- 1
    .Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
    while (.Object@binary_matrix[i, j] > 1) {
      .Object@binary_matrix[i, j] <- 0
      j <- j + 1
      .Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
    }
    .Object@sum_ones[i] <- sum(.Object@binary_matrix[i,])
  }
  
  # Initialize the vector m
  .Object@m <- array(0, c(.Object@N))
  for (j in 2:.Object@N) {
    if (.Object@sum_ones[j] < 1 / .Object@alpha) {
      if((1 + .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
      else .Object@m[j] <- -1
    }
    else if (.Object@sum_ones[j] == 1 / .Object@alpha) .Object@m[j] <- 0
    else {
      if((.Object@nb_var + 1 - .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
      else .Object@m[j] <- -1
    }
  }
  .Object
})

setMethod(f = "show", signature = "AEP.frame", definition = function(object) {
  cat("\n@nb_obs = number of observations = ", object@nb_obs, "\n\n")
  cat("@nb_var = number of variables = ", object@nb_var, "\n\n")
  cat("@data = [not displayed]\n\n")
  #print(object@data)
  cat("@N = pow(2, nb_var) = ", object@N, "\n\n")
  cat("@binary_matrix = \n")
  print(object@binary_matrix)
  cat("\n@sum_ones = sum of each row of the binary matrix = ", object@sum_ones, "\n\n")  
  cat("@alpha = ", object@alpha, "\n\n")
  cat("@iteration = ", object@iteration, "\n\n")
  cat("@m =", object@m, "\n\n") 
  cat("@bandwidth =", object@bandwidth, "\n")
})

setMethod(f = "summary", signature = "AEP.frame", definition = function(object) {  
  cat(object@nb_obs, "observations of", object@nb_var, "variables \n")  
  cat("  alpha: ", object@alpha, "\n")
  cat("  #iteration: ", object@iteration, "\n")  
  cat("  variation in variable (bandwidth):", object@bandwidth, "\n")
})

is.aep.frame <- function(frame) class(frame)=="AEP.frame"

as.aep.frame <- function(data, alpha=2/(dim(data)[2]+1), iteration=1, bandwidth=rep(0, dim(data)[2])) {
  if (missing(data)) stop("syntax error")  
  stopifnot(is.data.frame(data), is.numeric(alpha), is.numeric(iteration), is.numeric(bandwidth))    
  if (dim(data)[2] == 1) stop("'data' must have at least two columns" )
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")  
  if (iteration < 1) stop("the parameter 'iteration' must be >= 1")  
  if (length(bandwidth) != dim(data)[2]) stop("the parameter 'bandwidth' must be of the same length as the number of column of 'data'")
  new("AEP.frame", data, alpha, iteration, bandwidth)
}

VHmeasure <-function(frame, hc) {  
  if (missing(frame) || missing(hc)) stop("syntax error")
  stopifnot(is.aep.frame(frame), is.hcsimplex(hc))
  if (length(hc@b) != frame@nb_var) stop("the parameter 'hc@b' must be of the same length as the number of variables in 'frame")      
  # variable to increase the execution speed
  U.upper <- t(apply(frame@data, 1, function(d, s) s - d, s=hc@b + hc@h))
  U.upper <- U.upper / frame@data.weighted_bandwidth
  U.upper[U.upper > 1] <- 1
  U.lower <- t(apply(frame@data, 1, function(d, s) s - d, s=hc@b))
  U.lower <- U.lower / frame@data.weighted_bandwidth  
  U.lower[U.lower < -1] <- -1  
  Kn <- vapply(1:nrow(U.upper), 
    function(i) {
      if (any(U.lower[i,] > 1, U.upper[i,] < -1)) 0 
      else prod(U.upper[i,] - 1/3 * U.upper[i,]^3 - U.lower[i,] + 1/3 * U.lower[i,]^3)
    }
    , NA_real_)
  (1/frame@nb_obs) * (3/4)^(frame@nb_var) * sum(Kn) 
}


resolveAEP <- function(frame, vsum) {
  if (missing(frame) || missing(vsum)) stop("syntax error")  
  stopifnot(is.aep.frame(frame), is.numeric(vsum))    
  prob <- 0
  nb_simplex <- 0
  spx <- new("AEP.HCSimplex", array(0, c(frame@nb_var)), vsum, 1, 1)
  ListOfSimplexes <- c(spx)
  while (prob < 1 & length(ListOfSimplexes) > 0 & nb_simplex < frame@iteration) {    
    spx <- ListOfSimplexes[[1]]    
    if (length(ListOfSimplexes) > 1) ListOfSimplexes <- ListOfSimplexes[2:length(ListOfSimplexes)]
    else ListOfSimplexes <- NULL
    hc_full <- new("AEP.HCSimplex", spx@b, spx@h, spx@s, spx@iteration)
    prob_full <- VHmeasure(frame, hc_full)		
    nb_simplex <- nb_simplex + 1
    if (prob_full > 0) {
      hc <- new("AEP.HCSimplex", spx@b, frame@alpha * spx@h, spx@s, spx@iteration)
      prob <- prob + hc@s * VHmeasure(frame, hc)			
      nb_simplex <- nb_simplex + 1
      new.it <- spx@iteration + 1
      # Check how many new simplex can be generated without exceeding the limit frame@iteration      
      num_of_simplex <- frame@iteration - (nb_simplex + length(ListOfSimplexes))      
      if (num_of_simplex > 0) {                   
        for (j in 2:min(frame@N, num_of_simplex)) {
          new.b <- spx@b + (frame@alpha * spx@h * frame@binary_matrix[j,])
          new.h <- (1 - (frame@alpha * frame@sum_ones[j])) * spx@h
          new.s <- spx@s * frame@m[j]				
          ListOfSimplexes <- c(ListOfSimplexes, new("AEP.HCSimplex", new.b, new.h, new.s, new.it))			
        }			  
      }			
    }    
  }      
  new("AEP.result", prob, vsum, nb_simplex)
} 


AEP.d <- function(data, vsum, alpha=2/(dim(data)[2]+1), iteration=1, bandwidth=rep(0, dim(data)[2])) {  
  if (missing(data) || missing(vsum)) stop("syntax error")  
  stopifnot(is.data.frame(data), is.numeric(vsum), is.numeric(alpha), is.numeric(iteration), is.numeric(bandwidth))  
  if (dim(data)[2] == 1) stop("'data' must have at least two columns" )
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")  
  if (iteration < 1) stop("the parameter 'iteration' must be >= 1")	
  if (length(bandwidth) != dim(data)[2]) stop("the parameter 'bandwidth' must be of the same length as the number of column of 'data'")
  frame <- as.aep.frame(data, alpha, iteration, bandwidth) 
  resolveAEP(frame, vsum)
}

AEP.f <- function(frame, vsum) {    
  if (missing(frame) || missing(vsum)) stop("syntax error")
  stopifnot(is.aep.frame(frame), is.numeric(vsum))  
  resolveAEP(frame, vsum)
}


# AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- 
# AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- 

setClass("AEP.result", representation(
  prob = "numeric",
  vsum = "numeric",
  nb_simplex = "numeric"
)
)
setMethod("initialize", "AEP.result", function(.Object, prob, vsum, nb_simplex) {
  .Object@prob <- prob
  .Object@vsum <- vsum
  .Object@nb_simplex <- nb_simplex  
  .Object
})
setMethod(f = "show", signature = "AEP.result", definition = function(object) {
  cat("P[ sum <= ", object@vsum ,"] = ", object@prob, " ( x ", object@nb_simplex , ")\n")
})

# AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- 
# AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- 

setClass("AEP.HCSimplex", representation(
  b = "array",
  h = "numeric",
  s = "numeric",
  iteration = "numeric"
  )
)
setMethod("initialize", "AEP.HCSimplex", function(.Object, or, lim, multiplier, it) {
  .Object@b <- or
  .Object@h <- lim
  .Object@s <- multiplier
  .Object@iteration <- it
  .Object
})

is.hcsimplex <- function(hc) class(hc)=="AEP.HCSimplex"

# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 
# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 
# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 

  par(ask=FALSE)
	par(mfrow=c(1,1))
	
  my.data <- read.table("analysis_traces.txt", header=TRUE, sep=" ")

  my.nb_obs <- dim(my.data)[1]
  my.nb_var <- dim(my.data)[2]
  
  my.nb_points <- 100
  my.nb_iterations <- 500  
  #my.bandwidth <- c(0,0,0,0,0)
  my.bandwidth <- c(0.02,0.03,0.04,0.05,0.01)
  my.epsilon <- 0.0001
  my.max.range = 2

  # mydata cannot contain any 0 !!!!!
  # otherwise the estimation is screwed because of the P[X_k  in (b_k, b_k + h]] in equation 1.3. The problem is the "(" with b_k not included.
  # try with 2 variables X_1 and X_2 and the following values:
  # (10, 0) and (0, 50)
  # For Q( (0,0) , 50 ) we get H(0,0) - H(0, 50) - H(50, 0) + H(50, 50) = 0 - 0.5 - 0.5 + 1 = 0 !
  
  my.data[my.data == 0] <- my.epsilon

	observed.exectime <- array(0, c(my.nb_obs))
	for (i in 1:my.nb_obs) observed.exectime[i] <- sum(my.data[i,])
  
	#plot( x = 1:my.nb_obs, y = observed.exectime, 
	#	main = "Observed execution time",
	#	xlab = "#run",
	#	ylab = "Observed Exec Time",    
	#	type = "l", 
	#	col = "red")

	ECDF <- ecdf(observed.exectime)
  
  observed.ecdf <- array(0, c(my.nb_points))
  
  
  index <- 1  
  my.steps <- seq(min(observed.exectime), my.max.range * max(observed.exectime), length.out = my.nb_points)  
  for (i in my.steps) {
    observed.ecdf[index] <- ECDF(i)
    index <- index + 1
  }
  
	
  
  print("now processing V-AEP")
	
	estimated.ecdf <- array(0, c(my.nb_points))
  
	my.aep.frame <- as.aep.frame(data=my.data, iteration=my.nb_iterations, bandwidth=my.bandwidth)
  
	index <- 1
	for (i in my.steps) {
    result <- AEP.f(frame=my.aep.frame, vsum=i)
	  #result <- AEP.d(data=my.data, iteration=my.nb_iterations, bandwidth=my.bandwidth, vsum=i)
    print(result)
		estimated.ecdf[index] <- result@prob
		index <- index + 1
	}
	
  for (i in seq(length(estimated.ecdf), 1, -1)) {
    if (estimated.ecdf[i] > 1) estimated.ecdf[i] <- 1
    if (i < length(estimated.ecdf) & estimated.ecdf[i] > estimated.ecdf[i+1]) estimated.ecdf[i] <- estimated.ecdf[i+1]
  }

observed.Recdf = 1 - observed.ecdf
estimated.Recdf = 1 - estimated.ecdf


find_quantile <- function(x, y, q, precision) {
  f <- approxfun(x, y, 
                 method = "linear", 
                 yleft=1, 
                 yright=0, 
                 rule = 2, 
                 f = 0, 
                 ties = mean
  )  
  q <- 1 - q
  v_left <- min(x)
  v_right <- max(x)
  v <- v_left + (v_right - v_left) / 2
  while(abs(f(v) - q) > precision)  {
    if (f(v) - q < 0) v_right <- v
    else v_left <- v
    x <- f(v)
    v <- v_left + (v_right - v_left) / 2  
  }
  v
}

plot(x=my.steps, y=observed.Recdf,
     main = "Observed 1-ECDF vs. 1-Estimated ECDF",    
     xlab = "Execution time",
     ylab = "P[ExecTime > Y]",
     ylim = c(0,1),   
     lwd = 3,
     pch = 1,
     type="l"
)

lines( x = my.steps, y = estimated.Recdf, 		
  type = "l", 
  lwd = 3,
  pch = 1,
	col = "red"
)
  

legend("topright",
  legend = c("Observed 1-ECDF", "Estimated 1-ECDF"), 
  col = c("black", "red"), 
  text.col = "black", 
  bg = "white",
  lty = c(1, 1), 
  lwd=c(3,3),
  title="Legend", 
  bty='o', 
  title.adj=0.50
)

abline.lwd <- 2

estimated.q95 <- find_quantile(x=my.steps, y=estimated.Recdf, q=0.95, precision=0.001)
abline(v=estimated.q95, col="blue", lwd=abline.lwd)
text(estimated.q95, 0.9, " q = 95%", lty=3, col = "gray40", adj = c(0, -.1))

estimated.q99 <- find_quantile(x=my.steps, y=estimated.Recdf, q=0.99, precision=0.001)
abline(v=estimated.q99, col="blue", lwd=abline.lwd)
text(estimated.q99, 0.8, " q = 99%", lty=3, col = "gray40", adj = c(0, -.1))

observed.q95 <- find_quantile(x=my.steps, y=observed.Recdf, q=0.95, precision=0.001)
abline(v=observed.q95, col="green", lwd=abline.lwd)
text(observed.q95, 0.5, " q = 95%", lty=3, col = "gray40", adj = c(0, -.1))

observed.q99 <- find_quantile(x=my.steps, y=observed.Recdf, q=0.99, precision=0.001)
abline(v=observed.q99, col="green", lwd=abline.lwd)
text(observed.q99, 0.4, " q = 99%", lty=3, col = "gray40", adj = c(0, -.1))
